All-electron theory of the coupling between laser-induced coherent phonons in bismuth 
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Using first principles, all-electron calculations and dynamical simulations we study the behavior of the A\g 
and Eg coherent phonons induced in Bi by intense laser pulses. We determine the potential landscapes in the 
laser heated material and show that they exhibit phonon-softening, phonon-phonon coupling, and anharmonici- 
ties. As a consequence the Eg mode modulates the A\g oscillations and higher harmonics of both modes appear, 
which explains recent isotropic reflectivity measurements. Our results offer a unified description of the different 
experimental observations performed so far on bismuth. 

PACS numbers: 78.70.-g,63.20.Kr,71.15.Pd,63.20.Ry 



Intense femtosecond laser pulses can induce a nonequilib- 
rium state which leads to sudden and dramatic changes in the 
potential energy surfaces of different solids yj. This prop- 
erty can be used to excite and manipulate coherent lattice vi- 
brations, as has been recently shown by several experiments 

2e Ai 3\ and simulations L^. In the last years, a variety 

of experiments on laser excitation of coherent phonons has 

1^ n 

been performed on bismuth |7, 8, 9, 10], which is a particu- 
larly interesting solid since its ground-state structure exhibits 
a Peierls distortion. From the different studies done so far, a 
number of fundamental aspects still remain unexplained, like 
the detection of higher harmonics t9J and the appearance of 
modes that are forbidden by symmetry in isotropic reflectivity 
measurements |8]. 

In this letter we perform ab initio calculations which show, 
for the first time, the existence of a coupling between laser- 
induced phonon modes of different symmetry. In addition, we 
demonstrate that the large amplitude atomic vibrations excited 
by femtosecond laser pulses are affected by the anharmonic 
part of the potential energy surface, which creates overtones. 

The structure of Bi can be derived from a simple cubic 
atomic packing in two steps. First a simple cubic lattice is 
deformed by elongating it along one of the body diagonals, 
which is indicated by a thin line in Fig.^ A Peierls instability 




FIG. 1: (Color online) Structure of Bi (see text). In the Aig phonon 
mode the atoms move in the direction of the thin line as indicated for 
two atoms (long arrows). In the Eg phonon modes the atoms move 
in the perpendicular plane (e.g., in the directions of the short arrows). 



then causes the atoms to be displaced along the same diago- 
nal, in opposite directions (Fig.^. The magnitude of the dis- 
placements is determined by the detailed interaction between 
the Bi nuclei and the electrons. When the electrons are heated 
with a laser, the equilibrium positions of the atoms change. If 
the duration of the laser pulse is short in comparison with the 
timescale of the atomic motion, the ions are lifted at their orig- 
inal equilibrium positions to the potential energy surface that 
is created by the hot electrons, from where they start to swing 
about their new equilibrium positions 0] . This motion of the 
atoms corresponds to the coherent excitation of the displacive 
Aig phonons (Fig. 0. Coherent Eg phonons, which corre- 
spond to the motion of the atoms in the plane perpendicular 
to the elongated body diagonal, are also excited by a laser, 
through Raman scattering L^. In Bi the existence of laser- 
excited coherent Eg phonons has been demonstrated through 
electro-optical measurements fTTl . 

Recent experiments have shown three exciting results: (/) 
Through a measurement of the geometrical structure factors 
of two x-ray diffraction peaks the most accurate value so far 
of the maximal displacement of the atoms due to the inter- 
action with a short laser pulse has been determined 0]. {ii) 
In isotropic reflectivity measurements apart from the softened 
Aig phonon mode (2.52 THz) a signal at lower frequency 
(1.61 THz) that can be attributed to the Eg phonons has been 
detected 1 8]. In addition there is also a peak at 3.44 THz iV^ . 
(Hi) In another study higher harmonics of the softened Aig 
phonon frequency have been detected |9]. 

It is the aim of this letter to provide a theoretical explanation 
of these experimental results. To achieve this we performed 
a dynamical simulation for the Eg and Aig degrees of free- 
dom. The potential energy surface on which the atoms moved 
was calculated with the all-electron full-potential linearized 
augmented plane wave (LAPW) computer program WIEN2k 
II3ll . which is based on density functional theory (DFT). 

For our ground state DFT calculations we used the exper- 
imental unit cell parameters of Bi (spacegroup R3m, a = 
4.5332 A, c = 11.7967 A). The Bi atoms occupied the 6c 
sites (0,0,z), where z gives the position of the atoms along the 
above-mentioned diagonal. We performed DFT calculations 
for z = 0.230, 0.231, . . . , 0.250 c, where z = 0.250 c repre- 
sents the structure in the absence of the atomic displacement 
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due to the Peierls instability and z = 0.234 c is the experi- 
mental equilibrium structure. Based on the total energies for 
these z we could perform a simulation of the coherent A^g 
phonons. The degrees of freedom corresponding to the Eg 
phonons (x and y) were taken into account by varying x = 
0.000, 0.001, . . . , 0.010 a and z = 0.234, 0.235, . . . , 0.250 c. 
The potential energy surface, on which the ions moved, was 
then fitted to the function 

4 2 

(1) 

where the last term describes a coupling between the motion 
in the z and the x (y) directions. The fact that it is proportional 
to (a;^ + j/^ ) can be easily understood from the symmetry of the 
crystal structure, since the coupling cannot depend on the sign 
of the displacement (see Fig. [0. Other important details of 
our DFT calculations were as follows. LAPW's with energies 
up to 12.0 Ry were included in the basis. Inside the atomic 
spheres (with radii of 1.376 A) additional 5c/, 6s, 6p, and 6d 
local orbitals were used. Spin-orbit coupling was treated in 
a second variational procedure, where the scalar relativistic 
eigenstates up to 3.0 Ry and local 6/91/2 orbitals were used as a 
basis for the relativistic calculation. The entire Brillouin zone 
was sampled with 512 k points using temperature smearing 
(Te = 1 mRy). 

At elevated temperatures we assumed that there was no ex- 
change of heat between the electrons and the ions. As a con- 
sequence the entropy of the electrons was a constant of mo- 
tion 1 14]. We further assumed that the electrons were per- 
fectly thermalized at all times due to electron-electron in- 
teractions. Therefore the occupation numbers of the Kohn- 
Sham states were always given by a Fermi-Dirac distribution. 
Recently, another approach has been proposed 
where the occupation numbers for the electrons and holes of 
laser-excited Te fisi [Till and Bi |10] are modeled with two 
temperature distributions using different chemical potentials 
for the electrons and the holes while keeping the number of 
electron-hole pairs constant. Although this so-called two- 
chemical-potential model might work well for semiconduc- 
tors we do not believe that it is appropriate for Bi at high laser 
fluences as we discuss below. 

To compare the predictions of the two-chemical-potential 
model and the model where the constant entropy of the elec- 
trons was the only constraint we performed additional calcu- 
lations assuming two chemical potentials. In agreement with 
ilQi1 we chose the electron and hole temperatures — ^ 
0.5 eV at z = 0.234 c. No heat was assumed to be exchanged 
between the electrons or the holes and the ions. Therefore the 
entropy of the electrons and holes were constants of motion. 
In this respect we did not follow the method of 
where the temperatures of the electrons and the holes have 
been kept constant, while the atoms moved on the total en- 
ergy surface, because it yields incorrect forces 

The potential energy surface £'tot(7e) on which the ions 
moved at elevated electronic temperatures was calculated 



from 

Etot (Te) = Etot (gs) + A£;band , (2) 

where £'tot(gs) was the self-consistent total energy of the 
electronic ground state and Ai?band — £'band(7e) — 
£'band(gs). This approach is based on the interpretation of 
the Kohn-Sham energies as single-electron excitation ener- 
gies. In standard temperature-dependent DFT the electronic 
occupation numbers are incorporated in the self-consistent cy- 
cle to take into account possible screening effects. We also 
performed such standard temperature-dependent DFT calcu- 
lations and found that the differences between the predictions 
of both approaches were very small (the main effect was that 
the electronic entropies in both approaches had to be scaled). 
For this reason and also because it is not clear that the self- 
consistent approach always leads to better results we used 
the computationally less time-consuming non-self-consistent 
treatment of Eq. Our results for different laser-induced 
initial electronic temperatures are summarized in tableU 

In the electronic ground state the total energy i^tot was 
minimized for z — 0.2346 c. At this parameter the Aig 
and Eg phonon frequencies were 2.89 and 1.94 THz, respec- 
tively, in reasonable agreement with experiment [z = 0.2341 
c, i/(ylig) = 2.92, and iy{Eg) = 2.16 THz] fl and with ear- 
Uer calculations [z = 0.2336 c, viAig) = 2.93 THz] \T^. 
The differences with the experimental results are probably due 
to the local density approximation |18], which we used for 
the exchange and correlation energy. The differences with the 
earlier calculations ] 10] are probably due to the Bi pseudopo- 
tential used in those calculations. The present all-electron cal- 
culations do not rely on a pseudopotential and should be more 
accurate. 

In the excited electronic states (table|l} the potential energy 
surfaces are flatter than the ground state potential energy sur- 
face (Fig. As a consequence the phonon frequencies are 
lower For large amplitude coherent phonons anharmonicity 
may further lower the phonon frequencies As mentioned 
above, in ]7] the maximal laser-induced displacement Az has 
been accurately determined: Under the influence of a laser 
pulse that leads to a softening of the Aig frequency from 2.92 
to 2.12 THz, Az = 0.0085 ± 0.0021 c. For the electronic 
entropy Se — 0.428 /cs/atom we reproduced this softening 
of the Aig phonon frequency. According to our calculations 
16% of the decrease of the Aig phonon frequency was due to 
the anharmonicity of the potential (The harmonic frequency 
was 2.24 THz). Our value for the amplitude of the Aig mo- 
tion was Az — 0.0075 c, in agreement with the experimental 
value ]?]. 

After the initial softening of the Aig phonons, the fre- 
quency returns to its original value within roughly 10 ps IsV 
We discuss the origin of this frequency hardening by com- 
paring our results and the experimental results of ]8], where 
the Aig phonon frequency has been resolved in time. For the 
electronic entropy Se = 0.300 fcs/atom {T^ = 17.6 mRy at 
z = 0.234 c) our calculated Aig phonon frequency was 2.45 
THz, the initial Aig frequency obtained in After 0.3 ps 
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TABLE I: Coefficients for the potential energy surfaces of Eq. Q as a function of the electronic entropy Se- In Eq. Q Etot is given in mRy 
/ atom, X and y are given in units of the lattice parameter a, and z is given in units of c. The electronic temperature Te and the number of 
electron-hole pairs N^-h (expressed in % of the valence electrons), which are not constants of motion, are given for z = 0.234 c. 
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FIG. 2; (Color online) Potential energy surfaces for the ground state 
(Te = 1 mRy) and for excited electronic states (Te = 17.6 and 22 
mRy at z = 0.234 c). Vertical arrows show transitions of the atoms 
from the ground state to an excited energy surface. The Aig and Eg 
phonon modes are indicated. On the surface labeled Te — 22 mRy 
the curvature in the x direction becomes negative when the atoms 
are at their maximal displacement in the z direction. This negative 
curvature is indicated by a light (yellow online) arrow. 



the Aig frequency in |8] increased to 2.52 THz. Our results 
indicated that at most 0.03 THz of this increase can be ex- 
plained by the anharmonicity of the potential energy surface 
(The harmonic frequency was 2.48 THz), in agreement with 
experiments using two time-delayed pump laser pulses fioll 
that have shown that the frequency of the Aig phonon mode 
averaged over five periods is the same within 1% independent 
of the amplitude of the oscillations for frequencies as low as 
2.65 THz. Note however that anharmonicity plays a substan- 
tial role in the phonon softening in |7] as discussed above. 
Instead a decrease of the electronic entropy to Se = 0.260 
ks/ atom yielding an Aig frequency of 2.52 THz is probably 
responsible for the frequency increase observed in |8]. The 
accompanying decrease in the number of electron-hole pairs 
by 15% (calculated at z = 0.234 c) corresponds to a decay 
time of 1.8 ps, which agrees well with the experimentally de- 
termined electronic background decay time of 1.78 ± 0.08 ps, 
which is independent of the laser fluence fl^ . The loss of 
entropy of the electrons near the surface is most likely due to 



diffusion of hot electrons away from the surface region (the 
penetration depth of the laser c? « 17 nm) and may also be 
partly due to the exchange of heat of the electrons with the 
ions via electron-phonon coupling, two processes that were 
not explicitly taken into account in our calculations. 

At this point it is necessary to discuss the limitations of the 
two-chemical potential model proposed in fl3\. According to 
J2O] the electron-hole recombination time in Bi at 300 K is 
about 1 ps. In a highly excited electronic state this time is 
expected to be considerably shorter. So within less than four 
phonon periods the electrons and the holes obtain a common 
chemical potential. To check the validity of the two-chemical- 
potential model we applied it and found that the chemical po- 
tential of the holes was greater than the chemical potential 
of the electrons, which means that upon reaching a common 
chemical potential the number of electron-hole pairs would 
increase, which would then give rise to a substantial further 
lowering of the phonon frequencies (we estimated that the 
ions would even go over the baiTier at z ~ 0.25 c), in con- 
tradiction with experiment |8]. Therefore, we believe that the 
two-chemical potential model | l^l does not provide a realistic 
description of the electron dynamics in Bi. We can however 
not exclude that this model might be appropriate for semicon- 
ductors or for Bi at lower laser fluences. 

We now consider the coupling between the Aig and the 
Eg phonons. Figures |3ja) and|3jb) show the intensity of the 
Fourier transform of the z coordinate of the atoms. Higher 
harmonics of the main Aig peak at 2.52 THz, which have 
also been observed experimentally |9], and which are a con- 
sequence of the anharmonicity of the potential, are indicated. 
For the initial velocity of the atoms in the x direction (coher- 
ent Eg phonons) we chose a value that gave a peak to peak 
amplitude Ax = O.4A2;, in agreement with 12 ill . In Fig. 
|3jb) it is clear from the peaks that are not higher harmonics 
of the main frequency that there is a considerable coupling 
between the Aig and Eg phonons. In Eq. Q the Aig phonon 
mode couples to + y^, a signal with double the Eg fre- 
quency, 2v{Eg). Accordingly an analysis of the x coordinate 
of the atoms showed that the frequency of the highest peak 
induced by the phonon-phonon coupling, labeled 2Eg in Fig. 
|3b), (3.32 THz) equaled 2v{E„). Experimentally this peak 
has been observed at 3.44 THz I12ll . So in our calculations 
v{Eg) was slightly lower than in the experiment 1 12], which 
is consistent with our ground state calculation. The frequency 
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FIG. 3: Fourier transforms (intensities) of the z coordinate of the 
atoms for two different simulations: a), b) Te = 16.2 and c), d) 
Te = 12.7 mRy at 2 = 0.234 c. All curves have been convoluted 
with a Gaussian with a full width at half maximum of 0.03 THz. The 
left panels a) and c) show that the height of the peak labeled "£5" 
depends very sensitively on the laser fluence, which was twice lower 
in c) than in a). Experimentally this has been shown in 



of the peak labeled ''Eg" in Fig.|3b) (1.72 THz) is given by 
2i'{Aig) — 2v{Eg). Experimentally this peak has been ob- 
served at 1.61 THz fs*]. Our value was slightly higher than in 
the experiment |8] because the calculated v{Eg) was slightly 
lower. Note that in this peak has been identified with the 
Eg mode. However, since there is no coupling term linear in 
X (y), and since the Eg mode produces modulations of the 
Aig oscillations, it is clear that no peak of the power spec- 
trum can have the frequency v{Eg), whereas the difference 
2v{Aig) — 2v{Eg) should be present. The fact that the "Eg" 
peak has an intensity of only I{"Eg") = 7 10^'^I{Ain) and 
that it has nevertheless been observed experimentally f8) in- 
dicates that the laser-induced amplitude of the Eg mode is 
probably larger than Aa; = 0.4Az, which we used in our sim- 
ulation. According to fllS the relative Raman cross sections 
of the Aig and Eg modes are indeed temperature dependent. 
The sensitivity of the results to the initial conditions is further 
underlined by a calculation with Aa; — 1.2Az, for which we 
found I{"Eg") = 3 10"^/(74ig). Finally, we note that the 
peak labeled A in Fig.|3b) [0.80 THz = 2iy{Eg) - iy{Aig)] 
is experimentally hard to observe because of broadened spec- 
tral structures near zero frequency |12]. Similarly the peak 
labeled F in Fig.^tb) [2iy{Eg) + iy{Aig)] probably has too 
low an intensity to be observed. 

Figures |3jc) and |3jd) were calculated for the case where 
the energy absorbed from the laser was twice as small as in 
Figs.EJa) andHb) (T^ = 12.7 mRy at z = 0.234 c). It shows 
that the intensity of the peak "Eg" became orders of magni- 
tude lower upon reducing the laser fluence by a factor of two. 



Therefore a strong laser pulse is needed to observe the Aig-Eg 
phonon coupling, which has also been found experimentally 
in|8|. 

Another interesting result that we found is that the cou- 
pling between the Aig and Eg modes became very strong for 
Se > 0.431fcB/atom (T^ = 22 mRy at z = 0.234 c): In 
this case the curvature of the potential energy surface in the x 
direction became negative every time z reached its maximum 
value (Fig. |2j. We found that the amplitude of the Eg mode 
increased exponentially even when the initial velocity in the x 
direction was very small. 

In conclusion, we presented a dynamical simulation of co- 
herent laser-induced Aig and Eg phonons in Bi on a potential 
energy surface that was obtained by all-electron, first prin- 
ciples calculations. Using these simulations we could explain 
the effects of anharmonicity and Aig-Eg phonon coupling that 
take place after the interaction with ultrashort, intense laser 
pulses. We found good agreement with all experiments per- 
formed so far QSBEllIlIll- 
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